This notebook presents a brief introduction to the mapping capabilities of R. First we will need to install a few packages: ggmap, sf, osmdata, and osrm which will allow us to retrieve and visualize spatial data. Warning: sf installation for MacOS or Linux requires a bit of extra work, please follow the instructions here.
install.packages("ggmap")
install.packages("osmdata")
install.packages("sf")
install.packages("osrm")
The package ggmap is essentially an extension of ggplot; it enables you to plot a map as the background of a ggplot graph. The package sf (simple features) provides a set of functions to manipulate spatial vector data. The package osmdata provides the functionality to download and use data from OpenStreetMap(OSM). The package OSRM (OpenStreetMap-Based Routing service) allows users to calculate travel times and routed distances (km) based on OSM network data.
library(tidyverse)
library(ggmap)
library(osmdata)
library(sf)
library(osrm)
Before we actually build a query, we need to know what are the features that are available in OSM. For every feature there will also be a corresponding set of tags. Let's take a glimpse of some of the possibilities for features and tags related to the "amenity" feature.
available_features() %>% glimpse()
"amenity" %>%
available_tags() %>%
glimpse()
There are plenty of categories to choose from, now let's build a query. The getbb( ) function uses an API provided by OpenStreetMap (OSM) and returns a "bounding box" associated to a given place. opq() stands for "overpass query" and will allow us to build our query. Criteria is filtered using the add_osm_feature( ) function. Our query (q) will look for cinemas in Christchurch. add_osm_feature() takes two arguments, a feature (key="amenity") and a tag (value="cinema).
q <- getbb("Christchurch")%>%
opq()%>%
add_osm_feature("amenity", "cinema")
What is stored inside of q? Our query is actually a list containing four character vectors including the geographhic coordinates of our bounding box, and our search criteria.
q %>% glimpse()
Let's begin to extract some data. The osmdata_( ) function sends the query to the server and, depending on the suffix sf/sp/xml, returns a simple feature, spatial or XML format. A simple feature is standard file that used to access and work with one or two-dimensional geometries (point, line, and polygon) used in Geographic Information Systems. Let's try to extract the data in the form of a simple feature, what is stored inside of cinema?
cinema <- q %>%
osmdata_sf()
cinema %>% glimpse()
The result is a list of different objects, including simple feature classes and we are going to use ggplot's geom_sf() function to visualize them. geom_sf() uses a geometry object as its unique aesthetic, so in order to avoid conflicts with the default plot specification we set inherit.aes to FALSE, please have a look at the documentation here. The two main plotting features of ggmap are the get_map() function, which downloads a background map from a specified map provider (i.e Google Maps, OSM, Stamen Maps), and the ggmap() function that plots the raster map provided by get_map(). If you want to get familiar with GIS terminology you can read about it here. Notice the maptype argument for ggmap; it allows to specify the map theme, and options available include: "satellite","terrain", "hybrid", etc. Once a map is downloaded and plotted, it is then possible to use the normal grammar of gglot to visually represent that data on the map.
chc_map <- "Christchurch" %>%
getbb() %>%
get_map(source = "stamen", maptype = "toner")
chc_map %>%
ggmap()+
geom_sf(data=cinema$osm_points,
inherit.aes =FALSE,
colour="green",
fill="blue",
alpha=.5,
size=5,
shape=25)+
labs(x="Longitude",y="Latitude")
There are other features available in OSM and they can be represented through other geometry types like lines. Let's create a new query and this time we will try to map cycle ways within Christchurch.
q_2 <- getbb("Christchurch")%>%
opq()%>%
add_osm_feature("cycleway", "lane")
cycleway <- q_2 %>% osmdata_sf()
chc_map %>%
ggmap()+
geom_sf(data=cycleway$osm_lines,
inherit.aes =FALSE,
colour="blue",
size=1) +
labs(x="Longitude",y="Latitude")
What if you wanted to find the shortest route between two points. The OSRM package that can take care of that. The osrmRoute gets the shortest path between two points, it sends a query and gets the geometry of a route. The main arguments are the origin and destination, and in our case we will be passing them as sf objects. We will use two cinema locations. The remaining arguments specify the level of resolution required for our geometry and the type of object to be returned, in our case we are interested in obtaining a simple feature class.
route <- osrmRoute(src=cinema$osm_points[1,],dst=cinema$osm_points[2,],
overview="full", returnclass="sf")
str(route)
Notice the contents of the route object, specifically, there are two variables containing the travel time in minutes and travel distance in kilometers. Let's plot our route. We will use ggmap as we did before.
chc_map %>% ggmap() + geom_sf(data=route, inherit.aes = FALSE, colour ="blue")
Think of isochrones as polygons that represent areas that fall within a certain travel time from a given location. osrmIsochrone is a function that enables users to build isochrone maps. It takes the location of a given point and a vector that specifies the number and magnitude of "travel time" breaks. We also specify the number of additional points that will eb used to calculate the isochrones.
iso <- osrmIsochrone(loc = cinema$osm_points[2,], returnclass="sf",
# breaks = seq(from = 0, to = 14, by = 2), res = 50)
str(iso)
Notice that our object iso contains a "center" variable which actually corresponds to travel times in minutes for different regions. We will use this variable within the aesthetics of geom_sf.
chc_map %>% ggmap() + geom_sf(data = iso, aes(fill = center), alpha = 0.8, inherit.aes = FALSE) +
geom_sf(data=cinema$osm_points[2,], color = "blue", size = 2, inherit.aes = FALSE)
It is time to work on more interesting stuff. The fires in the Amazon have recently drawn a lot of media attention. NASA currently runs the Fire Information Resource Management System (FIRMS) and they have really interesting data that we can use for our exercise. The csv file is provided, however you can download the data yourself from here. The file we will be working with contains data on active fires in South America for the past 7 days. It contains the geographic coordinates of the fires along with other indicators such as the Fire Radiative Power (for more info, see previous link).
fire_data <- "https://firms.modaps.eosdis.nasa.gov/data/active_fire/viirs/csv/VNP14IMGTDL_NRT_South_America_7d.csv" %>%
read_csv()
fire_data %>% glimpse()
That is a lot of information as it contains records for South America. We want to narrow down the assessment to the region of MATO GROSSO which is one of the areas that has been affected recently. We can get the bounding box for that region. getbb actually returns a matrix with the coordinates of our bounding box.
MG_bb <- "STATE OF MATO GROSSO" %>%
getbb()
MG_bb
Now we can filter our dataset using the coordinates from our bounding box, and we only include observations with a high level of confidence.
fire_data_MG <- fire_data %>%
filter(longitude > -61.63, longitude < -50.21) %>%
filter(latitude > -18.04, latitude < -7.35) %>%
filter(confidence == "high")
fire_data_MG %>% glimpse()
Now we convert our data frame into an sf (simple feature) object that aside the original columns will have an additional geometry component. Then, with st_crs, we assign a geographic coordinate system, in our case we use WGS84 (World Geodetic Survey 1984). If you want to learn more about coordinate systems and projections, please go here
fire_MG_sf <- fire_data_MG %>%
st_as_sf(coords=c("longitude","latitude"))
st_crs(fire_MG_sf)
st_crs(fire_MG_sf) <- 4326
st_crs(fire_MG_sf)
str(fire_MG_sf)
Note the additional geometry list-column which holds the coordinates of all of our points. Let's extract some more data. As you have seen previously, OSM uses lots of features and tags to describe different things within a given area. For this particular example, we are interested in land use, more specifically on farmland and protected areas. Let's build a few queries upon these ideas. Let's start by using "landuse" as a key and "farmland" as a tag.
q_fm <- MG_bb %>%
opq() %>%
add_osm_feature("landuse", "farmland")
MG_farmland <- q_fm %>% osmdata_sf()
# MG_farmland %>% glimpse()
For the next query we use "boundary" as a key and "protected area" as a tag.
q_pa <- MG_bb %>%
opq() %>%
add_osm_feature("boundary","protected_area")
MG_protected <- q_pa %>% osmdata_sf()
# MG_protected %>% glimpse()
Now we are ready to start making some plots. First we extract a background "satellite" map for our region of interest. Let's take a look at both, farm and protected areas.
MG_bg_map <- MG_bb %>%
get_map(maptype = "satellite")
MG_bg_map %>% ggmap()+
geom_sf(data=MG_farmland$osm_polygons,
inherit.aes =FALSE,
aes(fill="Farmland"),
alpha=.5,
size=0.1, show.legend = "polygon")+
geom_sf(data=MG_protected$osm_polygons,
inherit.aes =FALSE,
aes(fill="Protected Area"),
alpha=.5,
size=0.1, show.legend = "polygon")+
scale_fill_manual(values = c("Farmland" = "blue", "Protected Area" = "purple")) +
labs(x="",y="")
It looks like in some parts the farms and the protected areas are very close to each other, in some cases they almost overlap. What about the fires? Let's bring them to our plot. The function coord_sf allows us to define the coordinate system that will be used within our plotting functions, in this case we use it to ensure that all of the layers of our plot are based on the same coordinate system. We will assign a color aesthetic to our Fire Radiative Power (rfp) attribute.
MG_bg_map <- "STATE OF MATO GROSSO" %>%
getbb() %>%
get_map(maptype = "terrain")
ggmap(MG_bg_map) +
coord_sf(crs = st_crs(4326)) +
geom_sf(data=MG_farmland$osm_polygons,
inherit.aes =FALSE,
aes(fill="Farmland"),
alpha=.5,
size=0.1, show.legend = "polygon")+
geom_sf(data=MG_protected$osm_polygons,
inherit.aes =FALSE,
aes(fill="Protected Area"),
alpha=.5,
size=0.1, show.legend = "polygon")+
scale_fill_manual(values = c("Farmland" = "blue", "Protected Area" = "purple")) +
geom_sf(data=fire_MG_sf, aes(color = frp),inherit.aes = FALSE) +
scale_color_gradient(low="yellow", high="red")
This notebook a only a brief introduction to mapping in R. If you are more interested in these topics and want to dig deeper into geographic data analysis I suggest you check out the following Geocomputation with R . It was recently released, is open source and has lots of examples using New Zealand as a case of study.